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Abstract 

We  present  a  multiscale  model  for  numerical  simulation  of  dynam¬ 
ics  of  crystalline  solids.  The  method  couples  nonlinear  elastodynamics 
as  the  continuum  description  and  molecular  dynamics  as  another  com¬ 
ponent  at  the  atomic  scale.  The  governing  equations  on  the  macroscale 
are  solved  by  the  discontinuous  Galerkin  method,  which  is  built  up 
with  an  appropriate  local  curl-free  space  to  produce  coherent  displace¬ 
ment  field.  The  constitutive  data  are  based  on  the  underlying  atom¬ 
istic  model;  it  is  either  calibrated  prior  to  the  computation  or  obtained 
from  molecular  dynamics  as  the  computation  proceeds.  The  decision 
to  use  either  the  former  or  the  latter  is  made  locally  for  each  cell  based 
on  suitable  criteria. 
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1  Introduction 


The  conventional  computational  methods  for  solid  mechanics  have  been 
primarily  based  on  continuum  models,  where  one  uses  a  small  number  of 
variables,  such  as  the  strain  and  stress,  to  efficiently  describe  the  mechanical 
properties.  The  quality  of  these  continuum  models  depends  heavily  on  the 
constitutive  assumptions,  which  are  usually  obtained  from  experimental  ob¬ 
servations.  On  the  other  hand,  atomistic  models,  such  as  molecular  statistics 
and  molecular  dynamics  (MD),  which  account  for  detailed  crystal  structure 
and  atomic  conhgurations,  have  been  proven  to  be  a  useful  methodology. 
However,  to  model  a  macroscale  process,  such  atomistic  systems  are  too  large 
to  £t  in  a  realistic  computation.  Therefore,  it  is  computationally  advanta¬ 
geous  to  develop  a  hybrid  model  that  includes  both  components.  Multiscale 
methods  that  aim  to  bridge  different  scales  have  been  an  active  area  of  in¬ 
terest.  Recently  developed  methods  in  this  class  include,  for  instance,  the 
quasicontinuum  methods  [44,  32],  the  Macro  Atomistic  Ab  initio  Dynamics 
method  (MAAD)  [2],  the  bridging  scale  method  [47],  the  bridging  domain 
methods  [48],  the  heterogeneous  multiscale  method  (HMM)  [23,  36,  24],  and 
the  pseudo-spectral  method  [45].  The  advantage  of  such  concurrent  methods 
is  that  the  models  can  be  refined  as  the  computation  proceeds. 

This  paper  is  concerned  with  the  dynamics  of  elastic  solids.  We  will  use 
the  framework  of  HMM  [23,  36].  The  main  idea  of  HMM  is  to  formulate  both 
the  atomistic  and  continuum  models  in  terms  of  universal  conservation  laws, 
and  the  coupling  is  accomplished  by  ensemble  averaging.  The  implementa¬ 
tion  of  such  method  consists  of  three  components, 

1.  a  macro  solver  for  the  continuum  model, 

2.  a  micro  solver  to  equilibrate  the  atomistic  system  locally  to  the  appro¬ 
priate  ensemble, 

3.  an  averaging  procedure  to  obtain  data  that  are  needed  in  the  continuum 
model. 

The  main  advantage  of  this  formalism  is  that  one  is  able  to  use  the  macroscale 
equations  to  capture  the  elastic  waves,  and  the  microscale  models  are  em¬ 
ployed  as  a  supplemental  component  to  provide  accurate  constitutive  data, 
thereby  bypassing  any  empirical  model  for  the  equation  of  state. 

The  constitutive  relation  obtained  from  the  atomistic  models  are  typi¬ 
cally  nonlinear  and  temperature  dependent.  As  a  result,  shock  waves  may 
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develop,  which  give  rise  to  large  strain  and  velocity  gradients,  posing  a  great 
challenge  for  numerical  simulations.  In  addition  suitable  constraints  have 
to  be  imposed  to  form  a  coherent  displacement  held.  The  problem  is  fur¬ 
ther  complicated  by  the  presence  of  defects,  causing  large  local  deformation 
and  discontinuity  in  strain  and  displacement.  To  overcome  these  difficulties, 
we  will  make  use  of  the  discontinuous  Galerkin  (DG)  method,  a  hnite  el¬ 
ement  method  based  on  piecewise  polynomials  as  basis  functions  that  are 
completely  discontinuous  across  element  boundaries.  The  DG  method  will 
be  formulated  for  the  elastodynamics  problems  under  the  multiscale  setting. 

Stable  and  convergent  DG  methods  have  been  designed  for  linear  and 
nonlinear  partial  differential  equations  (PDEs)  including  hyperbolic  conser¬ 
vation  laws  [14,  16],  convection-diffusion  equations  [15],  elliptic  equations  [3], 
and  dispersion  wave  equations  [49].  Gomparing  with  traditional  Galerkin 
methodology  (e.g.  [28]),  the  advantage  of  the  DG  methods  includes  their 
hexibility  in  h-p  adaptivity  and  parallel  efficiency.  We  refer  to,  e.g.  [17]  for  a 
review  of  the  DG  methods.  For  our  purpose,  the  DG  methods  provide  their 
advantage  in  the  ability  of  resolving  sharp  wave  fronts,  and  they  can  ac¬ 
commodate  a  locally  curl-free  function  space  for  the  deformation  gradient  to 
provide  a  coherent  displacement  field,  following  the  local  structure  preserving 
DG  methodology  in  [12,  35].  We  remark  that,  in  the  recent  years,  there  have 
been  many  contributions  in  developing  DG  methods  for  solving  PDEs  in  solid 
mechanics.  In  most  cases,  attention  has  been  restricted  to  linear  elasticity 
problems.  For  example.  Riviere  et  al.  [41]  formulated  and  analyzed  a  DG 
method  for  linear  elasticity  based  on  a  generalization  of  the  nonsymmetric 
interior  penalty  DG  method  for  the  diffusion  equation.  Kaser  et  al.  [30,  22] 
proposed  a  DG  method  in  space  that  uses  the  arbitrary  high-order  deriva¬ 
tives  (ADER)  approach  for  the  time  integration  to  solve  the  linear  elastic 
wave  equation  in  heterogeneous  media  on  unstructured  meshes  for  both  two- 
and  three-dimension.  Huang  and  Gostanzo  [27,  19]  proposed  a  space-time 
DG  formulation  for  linear-elasticity  where  stress  discontinuities  were  consid¬ 
ered  through  jumps  in  the  material  properties,  see  also  [31].  Abedi  et  al.  [1] 
proposed  a  space-time  DG  method  for  linearized  elastodynamics  that  deliv¬ 
ered  exact  balance  of  linear  and  angular  momentum  over  every  space-time 
element. 

The  main  purpose  of  this  paper  is  to  present  a  multiscale  methodology 
for  transient  elastodynamics  problems  based  on  the  DG  methods  as  a  macro 
solver.  We  will  discuss  how  the  DG  methods  can  be  used  to  capture  large 
scale  elastic  held,  and  how  the  molecular  dynamics  can  be  used  at  the  atomic 
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scale  and  coupled  to  the  macroscale  DG  solver  to  provide  accurate  constitu¬ 
tive  data. 


2  Macroscopic  and  microscopic  models 

Our  computation  involves  models  at  both  macroscopic  and  microscopic 
scales:  elastodynamics  on  the  macroscopic  scale,  describing  the  evolution  of 
the  elastic  held,  and  molecular  dynamics  at  the  atomic  scale,  providing  the 
constitutive  data  based  on  detailed  atomic  interactions.  In  this  section,  we 
briehy  present  the  equations  underlying  these  models. 

On  the  continuum  scale,  the  governing  equations  for  the  elastodynamics 
are  a  set  of  PDEs.  To  begin  with,  we  hx  a  reference  coordinate,  denoted 
by  X,  and  after  deformation,  the  point  will  be  displaced  to  a  new  position, 
x-|-u(x,  t),  with  u  being  the  displacement.  Let  e  =  Vu  be  the  corresponding 
deformation  gradient.  Then  the  continuum  equations  take  the  form: 

f  =  0, 

)  p^v-V-P  =  0,  (2.1) 

I  Pofe  +  V-j  =  0. 

Here  v  and  e  are  the  velocity  and  specihc  energy  per  particle  respectively, 

Po  is  the  initial  density,  P  is  the  hrst  Piola-Kirchhoff  stress  tensor  and  j 
is  the  energy  hux.  The  hrst  equation  in  (2.1)  describes  the  time  evolution 
of  the  deformation;  the  second  and  third  equations  are  the  conservation  of 
momentum  and  energy  respectively.  In  continuum  mechanics,  for  instance 
[33] ,  these  equations  are  supplemented  by  the  empirical  constitutive  relations 
for  stress  and  energy  huxes.  The  purpose  of  the  present  paper  is  to  develop 
multiscale  strategies  that  bypass  these  empirical  constitutive  laws  when  their 
accuracy  is  in  doubt. 

At  the  atomic  scale,  the  motion  of  the  atoms  constituting  the  solids  is 
given  by  molecular  dynamics. 


f  qi  =  Pi/rrii,  .  . 

I  P*  = 

where  m*  denotes  the  mass  of  the  Ah  atom,  q*  and  pj  are  the  generalized 
coordinate  and  momentum  for  the  Ah  atom,  and  H(qi,q2,---  ,  qvr)  is  the 
interatomic  potential  that  models  the  interaction  among  the  atoms.  In  this 
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paper,  we  will  only  consider  pair  potential.  Namely, 

^  =  rij  =  |rij|,  Tij  =  qi-  q^. 

Many-body  potential  models,  such  as  the  embedded  atom  model  (EAM)  [21] 
and  Tersoff-Brenner  model  [46],  can  be  dealt  with  similarly.  The  system  (2.2) 
is  a  Hamiltonian  system  with  the  Hamiltonian 

^^(q.p)  =  ‘5“ (qi.qz.'"  .qw).p  =  (pi.pa,  - ■■  .pw). 

(2.3) 

The  basic  assumption  in  our  current  method  is  the  scale  separation. 
Namely  the  relaxation  time  for  the  microscopic  processes  is  much  shorter 
than  the  typical  time  scale  of  the  continuum.  In  this  case,  it  suffices  to  con¬ 
sider  the  local  equilibrium  distribution  of  the  atomistic  system.  In  the  case 
of  zero  temperature,  this  is  equivalent  to  the  Cauchy-Born  hypothesis,  which 
suggests  that  the  atomic  displacement  follow  the  macroscopic  deformation. 
As  a  result,  the  strain  energy  density  and  elastic  stress  can  be  computed  from 
the  atomistic  models.  At  hnite  temperature,  one  has  to  take  into  account 
the  thermal  fluctuation  and  compute  elastic  properties  based  on  statistical 
ensembles.  In  this  paper,  we  consider  two  statistical  equilibrium;  the  micro- 
canonical  distribution, 

pi(q,p)  =  (2.4) 

where  the  volume  and  the  energy  of  the  system  are  prescribed,  and  the 
canonical  distribution, 

P2(q,p)  =  13  =  {keT)-^  (2.5) 

where  the  volume  and  the  temperature  are  given.  The  constant  Z,  which 
is  the  partition  function,  normalizes  the  probability  density.  Although  these 
statistical  ensembles  are  in  principle  equivalent  in  the  inhnite  volume  limit, 
one  might  have  more  practical  convenience  than  the  other.  Once  the  equi¬ 
librium  distribution  is  available,  physical  observables  can  be  computed  via 
statistical  averaging.  More  specihcally,  let  w  be  any  observable  with  a  mi¬ 
croscopic  expression  to (q,  p).  Then 

(u;)  =  y'M;(q,p)p(q,p)dqdp.  (2.6) 
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In  a  molecular  dynamics  simulation,  we  can  replace  the  ensemble  average  by 
a  time  average, 

(w)  =  lira  ^  [  w{q{t),p{t))dt,  (2.7) 

T^oo  I  Jq 

provided  that  the  system  is  ergodic. 

The  remaining  step  in  the  atomic/continuum  coupling  is  to  dehne  the  cor¬ 
responding  macroscopic  quantities  at  the  atomic  scale.  Since  the  continuum 
equations  (2.1)  are  expressed  in  the  Lagrangian  coordinate,  we  first  define 
the  reference  coordinate  for  the  atoms  as  the  equilibrium  positions,  denoted 
by  Qi-  Next  we  will  mainly  focus  on  the  calculation  of  the  stress  tensor.  The 
hrst  approach  to  dehne  stress  from  atomic  scale  is  from  thermodynamics  [7]. 
For  the  atomistic  system,  the  free  energy  is  dehned  as, 

T  =  —ksT  la  [  e~^^dqdp. 


Implicit  in  the  integral  is  the  dependence  on  the  deformation  gradient  e, 
which  changes  the  volume  and  shape  of  the  entire  system.  Here  periodic 
boundary  condition  is  assumed.  From  the  second  law  of  thermodynamics, 
we  have, 

n  de 

hxing  the  temperature.  Here  is  the  volume  of  the  system  in  the  reference 
coordinate.  This  calculation  yields. 


P 


®  Qi 


Qij 


Qi  Q  j  • 


(2.8) 


A  more  natural  approach  to  derive  the  microscopic  expressions  is  based  on 
conservation  laws  at  the  atomic  scale  [29].  More  specihcally,  we  dehne  the 
local  momentum  and  energy 


f  q(x,t) 
I  e(x,t) 


^qi(t)(5(x-  Qi), 
i 


Qi 


(2.9) 


As  a  result,  the  total  momentum  and  energy  in  a  control  volume  can  be 
easily  computed.  Furthermore,  one  can  derive  conservation  laws  for  the  local 
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where, 


=  ®  (<5*  -  Qj) 

X  [  6{:s.- {Qj  +  X{Qi  -  Qj)))dX, 

Y 

^  4YI  +  PiMi)  •  -  ^j)  ®  iQi  -  Qj) 

i^j 

X  f  5(x- {Qj  +  X{Qi  -  Qj)))dX. 

Jo 


(2.11) 


Averaging  (2.11)  in  space,  one  also  arrives  at  (2.8). 

Equations  (2.1)  and  (2.10)  have  been  the  starting  point  of  HMM:  the 
continuum  and  atomistic  models  can  be  coupled  at  the  level  of  conservation 
laws.  Such  observation  is  also  reminiscent  of  the  concept  of  local  equilibrium 
in  the  kinetic  theory  for  gas  dynamics,  where  the  Maxwellian  distribution  can 
be  used  to  provide  the  equation  of  state,  and  therefore  close  the  continuum 
equations. 


3  Numerical  methodology 

3.1  Macroscale  solver:  the  DG  method 

In  this  paper,  the  discontinuous  Galerkin  (DG)  and  the  local  discontin¬ 
uous  Galerkin  (LDG)  methods  are  used  as  our  macroscale  solver.  We  will 
briefly  review  the  DG  method  for  solving  hyperbolic  problems  and  the  LDG 
method  for  solving  parabolic  problems.  For  more  details,  we  refer  the  reader 
to  the  series  of  papers  of  Gockburn,  Shu  and  their  coworkers  [14,  13,  11,  16], 
the  lecture  notes  [10]  and  the  review  paper  [17]. 

To  briefly  illustrate  the  ideas,  consider  a  one-dimensional  conservation 
law  on  a  given  interval  I  =  [a,b]'- 

ut  +  f{u)x  =  0. 
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We  divide  interval  I  into  N  cells  as  follows, 


a  =  xi  <  X3  <  . . .  <  Xj^.i  =  b. 

2  2  2 

(3.1) 

We  denote 

1 

Ii  =  {Xi_i,Xi^),  Xi  =  -{Xi_i  +  Xi^) 

(3.2) 

and 

Axi  =  Xi,i  —  Xi  1 ,  h  =  max  Ax*. 

2  2  i 

(3.3) 

Next  we 

dehne  the  approximation  space  as 

(3.4) 

Here  P^{Ii)  denotes  the  set  of  all  polynomials  of  degree  at  most  k  in  the 
interval  !{.  For  simplicity,  we  will  use  the  notations  u  and  v  instead  of  Uh, 
Vh  to  denote  the  numerical  solutions  whenever  they  do  not  cause  confusion. 
We  can  choose  a  basis  of  P^{Ii)  as,  for  example,  {1,  },  where  the 

monomials  f  For  a  function  v  G  we  use  v~  i  and  i  to  refer  to 

the  left  and  right  limit  of  v  at  respectively,  at  the  interface  where  v  is 
discontinuous. 

The  formulation  of  the  DG  method  is  as  follows:  hud  u{.,t)  G  such 
that  for  any  test  functions  v  E  V/^, 


f{u)v^dx  +  f{u).  iv  . 

2  2 


0. 


(3.5) 


The  single- valued  flux  f-j_i  should  be  taken  as  a  monotone  flux  depending 
on  both  u~_^i  and  (exact  or  approximate  Riemann  solvers  in  the  system 

case);  see,  for  example,  [34].  We  will  choose  the  numerical  fluxes  to  be  the 
Lax-Friedrichs  flux 


i+4 


(3.6) 


where  a  =  max  |/'(m)|. 

Time  discretization  is  done  by  the  nonlinearly  stable  high  order  TVD 
Runge-Kutta  methods  developed  in  [43].  In  particular,  the  second-  and 
third-order  TVD  Runge-Kutta  methods  are  used  to  match  the  corresponding 
spatial  accuracy  in  our  paper.  For  solving  the  method  of  the  lines  ODE 


Ut  =  L{u), 


(3.7) 


where  L{u)  can  be  any  spatial  discretization  of  u,  the  second-order  TVD 
Runge-Kutta  method  is  given  by 

=  u^  +  /\tL{u^)  (3.8) 

and  the  third-order  TVD  Runge-Kutta  method  is  given  by 

=  u^  +  /\tL{u^) 

^(2)  _  1^(1)  ^  (3_g^ 

4  4  4 

3  3  3  ^  ^ 

Since  the  DG  method  may  have  oscillations  when  the  solutions  contain 
discontinuities,  nonlinear  total  variation  bounded  (TVB)  limiters  are  often 
used.  The  limiter  we  use  in  this  paper  is  briefly  described  below.  For  more 
details  we  refer  the  readers  to  Shu  [42]  and  Cockburn  et  ah  [13,  16].  Denote 

u~^i=Ui  +  Ui,  u^_i=Ui-Ui,  (3.10) 


where  Ui  is  the  cell  average.  Ui  and  Ui  are  modihed  by 


u 


(mod)  /  ~  -  _  _  -  \  -(mod)  /  ~  _  _  _  -  \ 

A  )  =  m  [Ui,  UiJ^i  -  Ui,  Ui  -  Ui_i) ,  u]  '  =  m  [Ui,  -  Ui,  Ui  -  Ui_i)  , 

(3.11) 

where  the  minmod  function  m  is  given  by 


f  s-mini<j<„|aj|  if  sign(ai)=...=sign(a„)=s,  .  . 

0  otherwise. 


or  by  the  TVB  modihed  minmod  function  [42] 


?fi(ai,a2,  ...,a„) 


ai  if  |ai|  <  Mh‘^ 

m(ai,  02, ...,  a„)  otherwise. 


(3.13) 


where  M  >  0  is  a  constant.  The  choice  of  M  depends  on  the  solutions  of  the 
problem.  For  the  scalar  case,  it  is  possible  to  estimate  M  [M  is  related  to 
the  magnitude  of  the  second  derivatives  of  the  solution  at  smooth  extrema); 
however,  it  is  more  difficult  for  systems.  In  our  paper,  we  apply  the  limiter  to 
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every  component  of  the  system  with  a  suitable  M  (which  may  not  necessarily 
be  the  optimal  M). 

In  the  case  where  diffusion  terms  are  present,  e.g.  the  convection-diffusion 
problems 

ut  +  f{u)x  =  {a{u,  x)u^)^,  (3.14) 

where  a{x)  >  0,  the  idea  of  the  LDG  method  [15]  is  to  rewrite  (3.14)  into  a 
hrst  order  system 


Ut  —  {a{u,x)w)x  =  0,  w  —  Ux  =  0.  (3.15) 

We  can  then  formally  use  the  same  DG  method  for  the  convection  equation 
to  solve  (3.15),  resulting  in  the  following  scheme:  hnd  u{.,t),w{.,t)  G 
such  that 


aWi.iv.  1  -|-  aw,  1  =0 

*+2  2  ^-2 

(3.16) 

-  U,,  1Z~.+U,  IZ^  t  =  0 
*  2 

(3.17) 

for  all  test  functions  v,z  ^  In  [15],  criteria  are  given  for  these  fluxes  to 
guarantee  stability,  convergence,  and  a  suboptimal  error  estimate  of  order 
k  in  the  L2  norm  for  piecewise  polynomials  of  degree  k.  A  clever  choice  of 
fluxes  in  an  alternating  way 


or 


^i+7 


Ut+I  =  u 


y..  y+4=<. 

(3.18) 

y.,  *.+1=®-+^ 

(3.19) 

would  satisfy  these  criteria  and  give  a  scheme  of  order  k  +  1. 


3.2  Molecular  dynamics  simulation 

In  computing  the  elastic  stress  from  the  microscopic  model,  we  prepare 
the  atomistic  system  as  follows.  We  hrst  arrange  the  atoms  to  the  equilib¬ 
rium  position,  Qj,  and  then  apply  a  uniform  deformation,  =  {I  +  A)Qj. 
The  major  axes  of  the  simulation  box  are  also  deformed  accordingly.  Peri¬ 
odic  boundary  conditions  are  applied  with  respect  to  the  deformed  box  to 
maintain  the  deformation  gradient  throughout  the  computation.  The  initial 
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velocity  for  the  atoms  is  sampled  from  a  Gaussian  distribution  with  the  vari¬ 
ance  given  by  ksT.  In  order  to  maintain  the  system  at  a  given  temperature, 
we  use  the  standard  NosG Hoover  thermostat  [40,  26].  To  ensure  numerical 
stability,  the  time  step  is  chosen  so  that  5t  <  2/uJmax)  where  Umax  is  the 
largest  phonon  frequency.  During  the  simulation,  the  stress  is  sampled  ev¬ 
ery  20  time  steps,  and  then  averaged  at  the  end  of  the  simulation  according 
to  (2.7).  In  addition  the  Verlet  list  method  is  used  to  speed  up  the  force 
calculation.  These  standard  techniques  can  be  found  in  [25]. 

3.3  The  multiscale  method 

Having  described  the  methodology  for  solving  the  atomistic  and  contin¬ 
uum  equations,  we  are  now  in  a  position  to  describe  the  multiscale  method. 
The  coupling  strategy  is  quite  straightforward:  we  hrst  divide  the  computa¬ 
tional  domain  into  cells  over  which  the  macroscale  equations  (2.1)  are  dis¬ 
cretized.  In  computing  the  numerical  fluxes,  the  stress  and  heat  flux  are 
obtained  directly  from  the  atomistic  model,  either  from  a  molecular  dynam¬ 
ics  simulation,  or  from  a  simplihed  model,  also  calibrated  from  atomistic 
simulations  beforehand  (see  next  section).  This  procedure  can  be  demon¬ 
strated  from  Fig.  3.1.  Comparing  to  the  conventional  numerical  procedure, 
we  have  used  atomistic  model  as  a  supplemental  component  on  each  cell  to 
supply  the  constitutive  data. 

4  One-dimensional  shock  propagation 

In  our  hrst  example  the  system  is  described  at  the  atomistic  level  by  a 
Lennard-Jones  potential: 

All  the  particles  are  assumed  to  have  mass  m.  In  the  molecular  dynamics 
simulation,  the  parameters  e  and  a  are  normalized  to  unit.  As  a  result, 
all  the  computational  results  are  expressed  in  terms  of  reduced  units,  rep¬ 
resented  by  m,  e  and  a.  The  atomistic  system  is  two-dimensional  and  it  is 
constrained  to  a  narrow  slab  so  that  the  macroscopic  dynamics  is  essentially 
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one- dimensional.  The  continuum  equations  then  reduce  to 


{dtSii  —  dxVi  =  0, 

podtvi  -  dxPii  =  0,  (4.2) 

podte  -  dx{PiiVi)  =  0. 

The  initial  condition  for  the  macroscale  quantities  is  set  up  as  follows:  For 
the  a;  <  0  half  plane,  we  impose  a  uniform  deformation  gradient  eu  =  —0.01, 
and  for  the  other  half  plane,  the  deformation  gradient  is  zero.  The  system 
starts  from  zero  velocity  and  the  temperature  T  =  0.1  on  the  left  half  plane, 
and  T  =  0.3  on  the  right  half  plane.  From  the  continuum  viewpoint,  this  is 
an  example  of  the  Riemann  problem. 

We  apply  the  multiscale  procedure  described  above  to  this  problem  with 
the  DG  method  as  the  macroscale  solver.  The  P^  and  P^  (second  and  third 
order  based  on  piecewise  linear  and  quadratic  polynomials)  DG  results  are 
both  presented.  The  numerical  results  are  shown  in  Fig.  4.1.  As  a  com¬ 
parison,  we  also  compute  the  solution  with  the  second-order  Lax-Friedrichs 
type  central  schemes,  described  in  Nessyahu  and  Tadmor  [39].  At  t  >  0,  one 
observes  two  shocks  separated  by  a  contact  discontinuity.  We  can  see  that 
the  DG  method  captures  the  discontinuities  very  well.  This  type  of  results 
have  also  been  conhrmed  by  a  full  atomistic  simulation  [36] . 


5  Two-dimensional  case 

In  our  next  simulation,  the  atomistic  model  is  a  three-dimensional  Lennard- 
Jones  solid  with  the  face-centered  cubic  (FGG)  structure.  We  consider  a  sys¬ 
tem  in  a  state  of  plain  strain  so  that  the  macro  scale  behavior  is  essentially 
two-dimensional . 

5.1  Atomistic- based  constitutive  models 

The  energy  flux  j  in  the  continuum  equations  (2.1)  can  be  written  as 

j  =  -(v'^P  +  q),  (5.1) 


with  q  being  the  heat  flux. 

In  principle  the  heat  flux  can  be  computed  from  an  atomistic  simula¬ 
tion  with  a  consistent  temperature  gradient  [36].  However  such  simulations 
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typically  take  much  longer  time  to  relax,  making  the  computation  rather 
expensive.  As  an  alternative,  we  model  the  heat  flux  by  the  Fourier  law, 

q  =  kVT,  (5.2) 

which  has  been  confirmed  via  numerous  numerical  studies  [8] .  The  heat  con¬ 
ductivity  has  been  precomputed  from  the  atomistic  model.  In  the  following 
computations,  we  will  consider  both  the  case  with  the  heat  flux  as  well  as 
the  one  without  the  heat  flux. 

In  the  case  of  small  deformation  and  low  temperature,  a  linear  approxi¬ 
mation  can  be  made  to  provide  an  explicit  stress  strain  relation: 

Pij  CijkiEf^i  (yTSij  i)  j)  k,l  1)2  (5.3) 

where  Cijki  are  elastic  parameters.  They  can  be  computed  directly  from  the 
atomistic  model  [4].  In  particular,  we  have 

{Pn  =  +  <^12^22  +  ctT, 

Pu  =  Cu{e,2+e2i),  (5^^) 

fhl  —  ^12: 

P22  =  €11622  +  Ci2^ii  +  aT. 

In  terms  of  the  reduced  units,  these  parameters  have  been  obtained  from  the 
atomistic  model, 

Cii  =  97.56753,  C12  =  55.56524,  (^44  =  55.55364, 

a  =  -9.91139,  K  =  0.1925353265. 

It  has  been  calibrated  that  under  the  condition  that, 

||e||2  <  0.004  and  T  <  0.2,  (5.5) 

the  error  of  this  linear  approximation  (5.4)  is  within  one  percent.  However 
this  model  cannot  be  directly  used  because  in  the  PDFs  the  temperature  is 
not  a  derived  quantity  and  needs  to  be  obtained  from  the  atomistic  models. 
For  this  purpose  we  make  another  approximation  for  the  temperatnre, 

PqC  =  E{e,  T)  =  Eq  +  -CijkiSkiSij  +  3T.  (5.6) 
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Here  the  first  term  is  the  cohesive  energy  of  the  crystalline,  Eq  =  —7.4591 
per  atom,  the  second  term  is  the  familiar  form  of  the  strain  energy  for  linear 
elasticity.  In  the  case  of  low  energy,  namely, 

\E-Eo\<l.2,  (5.7) 

the  modeling  error  for  this  approximation  (5.6)  is  also  within  one  percent. 

The  equations  (5.4)  and  (5.6)  provide  an  efficient  constitutive  model  in 
the  regime  of  small  strain  and  low  temperature.  It  can  be  shown  that  with 
these  constitutive  equations,  the  system  of  PDEs  is  hyperbolic  (see  Appendix 
A).  With  the  heat  flux  q  present,  the  system  is  of  parabolic  type. 

5.2  Implementation  of  the  constitutive  models 

While  the  equilibrium  molecular  dynamics  in  principle  offers  the  correct 
constitutive  equation,  such  procedure  is  typically  computationally  intensive, 
making  the  overall  computation  rather  expensive.  On  the  other  hand,  the 
simplihed  models  in  Sec.  5.1  are  efficient,  but  the  accuracy  can  only  be  guar¬ 
anteed  under  small  deformation  and  low  temperature.  Therefore  it  is  com¬ 
putationally  feasible  to  include  both  these  models  in  the  simulation:  we  use 
the  constitutive  relation  (5.4)  if  the  criterion  (5.5)  is  met,  and  (5.6)  if  the 
condition  (5.7)  is  satished;  Otherwise  the  data  are  obtained  directly  from  a 
molecular  dynamics  simulation.  This  reminds  us  to  apply  the  Domain  De¬ 
composition  Method  (DDM)  [9].  The  main  idea  of  the  DDM  is  to  solve  the 
relatively  inexpensive  macroscopic  model  in  most  part  of  the  computational 
domain,  and  only  solve  the  micro  problems  in  the  sub-domains  where  the 
macroscopic  model  is  not  valid. 

More  specihcally,  our  computational  domain  can  be  decomposed  to  four 
types  of  sub-domains.  These  sub-domains  are  distinguished  by  whether  the 
conditions  (5.5)  and/or  (5.7)  are  satisfied.  For  example,  the  first  type  of 
sub-domains  can  be  the  one  satisfying  both  (5.5)  and  (5.7).  In  this  case, 
a  computationally  inexpensive  macro  model  is  sufficient,  and  we  no  longer 
need  to  do  MD.  The  procedure  of  using  the  DDM  is:  at  each  time  step  and 
in  each  cell,  if  the  energy  is  low,  namely,  if  the  condition  (5.7)  is  satished, 
we  will  compute  the  temperature  from  (5.6)  instead  of  the  MD  procedure. 
After  the  temperature  is  obtained,  we  compute  the  stress  and  energy  hux 
from  the  atomistic  model  via  a  canonical  ensemble.  For  small  deformation 
and  low  temperature,  i.e.  when  (5.5)  holds,  we  use  (5.4)  to  bypass  the  MD 
procedure. 
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In  practice,  we  use  parallel  computing  to  speed  up  our  program.  At  the 
macro  computational  domain,  for  each  time  step,  we  hrst  need  to  get  the 
information  of  temperature.  We  label  the  L  cells  to  be  those  which  satisfy 
the  condition  (5.7),  and  the  H  cells  to  be  those  in  which  this  condition  is 
not  satished.  Then  we  use  the  macro  relation  (5.6)  to  get  the  temperature 
in  the  L  cells.  This  needs  only  a  negligible  cost  comparing  to  MD.  For  the 
H  cells,  we  do  need  to  use  MD  with  the  canonical  ensemble  (2.5),  which  are 
computationally  expensive.  To  efficiently  use  all  the  processors,  we  collect 
all  the  H  cells  and  distribute  them  to  available  processors  evenly.  Next,  we 
need  to  get  the  information  of  the  stress  tensor.  The  procedure  is  the  same 
as  above,  except  that  now  the  condition  (5.5)  and  the  relation  (5.4)  are  used 
to  distinguish  different  types  of  cells.  This  procedure  is  repeated  for  every 
macro  time  step. 

5.3  Curl-free  discontinuous  Galerkin  method 

In  the  hnite  volume  representation  of  the  elastic  held,  it  is  important 
to  ensure  the  coherent  structure  at  the  cell  interface.  This  issue  has  been 
recognized  on  [1] .  In  fact  the  macro  scale  PDFs  are  equipped  with  a  natural 
constraint: 

V  X  e  =  0, 

which  leads  to  two  pairs  of  curl-free  variables,  (£11,612)  and  (£21, £22)-  In 
order  to  fulhll  the  constraints,  we  choose  a  curl-free  basis  for  our  DG  method 
instead  of  the  standard  polynomial  basis.  The  original  idea  of  using  curl-free 
basis  is  from  [12],  in  which  the  locally  divergence- free  DG  method  is  designed 
to  solve  the  Maxwell  equations,  and  from  [35],  in  which  the  locally  curl- free 
DG  method  is  designed  to  solve  the  Hamilton-Jacobi  equations. 

In  the  standard  RKDG  method,  we  seek  the  solution  in  the  hnite  dimen¬ 
sional  polynomial  space 

=  {v  =  {Vi,V2,V3,V4,V5,Vg,Vy)  1  v\k  G  G  IC}  (5.8) 

where  K  is  the  element  and  /C  is  the  collection  of  the  elements,  P^’^(A')  = 
{P^{K)y,  and  P^{K)  denotes  the  space  of  polynomials  in  K  of  degree  at 
most  k.  Now  we  are  looking  for  the  solution  space  whose  bases  satisfy  the 
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curl-free  condition,  i.e. 


■yfc  _  J  a;-  c  •  ^'*^1  _  9v2  dv?,  _  dv^  I 

^  h  \  ^  h  '  dy  dx  ^  dy  dx  j 

e  !v=(t,3,t.4)  ^  =  1^1 

©  jv  =  (^;5,U6,^'7)  :  V  G 

=  ©  vf 

That  is,  the  vectors  (^1,^2)  and  (^3,^4)  are  curl-free  polynomial  vectors. 
Notice  that  the  dimension  of  the  space  is  {k  +  1)(A;  +  4)/2,  only  about 
half  as  the  dimension  of  which  is  (A;  +  l)(A;  +  2).  Thus  the  total  dimension 
of  the  curl  free  space  is  (/c^  +  k)  less  than  the  stand  piecewise  polynomial 
space  V^’^,  and  we  can  save  a  lot  of  computational  cost  by  using  the  curl 
free  space. 

It  is  very  easy  to  obtain  the  local  bases  for  V^’^.  We  can  take  the  gradient 
of  the  standard  bases  of  since  we  know  V  x  V/  =  0  for  any  scalar 

function  /.  For  example,  if  K  is  the  rectangle,  with  the  center  {xi,yj)  and 
width  Axi,  Ai/j,  if  we  denote 

x-Xi  y-Vj 

Axi  ’  ^  Ayj 

then  one  set  of  bases  of  would  be,  when  k  =  1, 

For  /c  =  2,  we  need  to  add 

f  e\  (  2A|/,er/  \  (  Ay,y^  W  0  \ 

0  y  ’  Axi^‘^  )  ’  2Axiiy  /  ’  V  / 

And  for  A;  =  3,  we  need  to  add 

(  e  \  (  \  (  \  (  A9,,,=  W  0  \ 

f  0  t  J  ’  t  3Ax,^/7^  )  '  \  rj^  )  ' 


(5.10) 

(5.11) 

(5.12) 
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5.4  Numerical  results 


5.4.1  Comparison  of  RKDG  methods  using  curl-free  and  stan¬ 
dard  P^  bases 

In  the  following  example,  we  consider  a  system  under  an  external  force 
\i{x,y,t).  The  macroscale  equations  take  the  form  of: 

r  =  0’ 

{  dofv-V-P  =  hi{x,y,t),  (5.13) 

I  dof  e  -  V  ■  (v^P)  =  h2{x,y,t). 

We  choose  the  source  terms  as 

hi{x,y,t)  =  0.01(po-C'ii-C'i2-2C44)cos(t  +  a;  +  |/), 
h2{x,y,t)  =  0.01^(po  -  C*!!  -  <^12  -  2C44)  sin(2(t  +  a;  +  I/)). 

The  initial  conditions  are  given  by, 

{Eij  =  0.01  sin(a;  +  I/) 

Vi  =  0.01  sin(a;  +  y)  i,j  =  l,2  (5.14) 

Pqc  =  Eq  +  0.01^(Cii  +  C12  +  2C44)  sin^(x  +  y) 

and  with  periodic  boundary  conditions  the  exact  solutions  are 

{Eij  =  0.01  sin(t  +  X  +  y) 

Vi  =  0.01  sin(t  +  a;  +  y)  i,j  =  l,2  (5.15) 

Pqc  =  Eq  +  0.01^(Cii  +  C12  +  2C44)  sin^(t  +  a:  +  p) 

We  use  this  set  of  exact  solutions  to  test  the  accuracy  and  efficiency  of 
the  curl-free  DG  method  comparing  to  the  standard  DG  method.  We  list 
the  errors  and  orders  of  accuracy  using  P^  and  P^  DG  methods  in  Tables 
5.1,  5.2  and  5.3.  The  solutions  are  run  to  T  =  27r.  We  can  observe  the 
optimal  (/c  -|-  l)-th  order  of  accuracy  both  for  the  standard  DG  method  and 
the  curl-free  DG  method.  The  magnitudes  of  the  errors  for  the  same  mesh 
are  comparable  for  the  two  DG  methods. 
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Table  5.1:  error  and  order  of  accuracy  for  the  stress  strain  e 

mesh  standard  local  curl-free 

error  order  L^error  order 


pi 

10  X  10 

2.12E-03 

- 

2.12E-03 

~ 

20  X  20 

5.28E-04 

2.01 

5.27E-04 

2.00 

40  X  40 

1.31E-04 

2.00 

1.31E-04 

2.00 

80  X  80 

3.32E-05 

1.98 

3.32E-05 

1.98 

160  X  160 

8.35E-06 

1.99 

8.35E-06 

1.99 

p2 

10  X  10 

2.35E-03 

- 

1.92E-03 

~ 

20  X  20 

3.25E-04 

2.86 

2.62E-04 

2.87 

40  X  40 

4.41E-05 

2.88 

3.39E-05 

2.95 

80  X  80 

5.71E-06 

2.95 

4.33E-06 

2.97 

160  X  160 

6.93E-07 

3.04 

5.47E-07 

2.98 

5.4.2  Thermal  expansion 

In  this  example,  we  study  the  effect  of  thermal  expansion  due  to  the 
temperature  dependence  of  the  stress.  Initially  the  material  is  at  rest  with  a 
homogeneous  temperature  distribution  T  =  0.1.  Our  computational  domain 
is  a  square  [—1,1]  x  [—1,1].  We  then  increase  the  temperature  in  the  middle 
[—0.4,  0.4]  X  [—0.4,  0.4]  instantaneously  to  T  =  0.4.  This  results  in  a  thermal 
expansion  that  propagates  outwards.  We  solve  the  macro  equations  (2.1) 
by  the  DG  methods  described  above,  and  compare  the  results  with  those 
obtained  with  the  two-dimensional  central  scheme  [39]. 

We  present  results  in  the  following  two  cases:  one  is  without  the  heat 
flux,  the  other  is  with  the  heat  flux.  Without  the  heat  flux,  the  system 
is  hyperbolic.  Therefore  it  will  have  discontinuities  in  the  solution  under 
discontinuous  initial  conditions.  In  the  presence  of  the  heat  flux,  the  system 
is  parabolic.  Therefore,  the  solution  becomes  smooth  at  t  >  0  even  if  the 
initial  condition  is  discontinuous.  We  will  demonstrate  the  numerical  results 
for  the  temperature  distribution  as  well  as  the  velocity  held  below. 

The  simplified  constitutive  model 

We  hrst  run  the  simulation  using  the  linear  constitutive  model  (5.4)  and 
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Table  5.2:  error  and  order  of  accuracy  for  the  velocity 

mesh  standard  local  curl-free 

error  order  L^error  order 


pi 

10  X  10 

2.29E-02 

- 

2.28E-03 

~ 

20  X  20 

7.88E-03 

1.54 

7.88E-03 

1.53 

40  X  40 

1.96E-03 

2.01 

1.96E-03 

2.01 

80  X  80 

4.70E-04 

2.06 

4.69E-04 

2.06 

160  X  160 

1.14E-04 

2.04 

1.14E-04 

2.04 

p2 

10  X  10 

2.56E-03 

- 

2.49E-03 

~ 

20  X  20 

3.03E-04 

3.08 

2.95E-04 

3.08 

40  X  40 

3.74E-05 

3.02 

3.68E-05 

3.00 

80  X  80 

4.65E-06 

3.01 

4.63E-06 

2.99 

160  X  160 

5.77E-07 

3.01 

5.82E-07 

2.99 

(5.6)  everywhere. 

Without  the  heat  flux;  The  results  of  the  temperature  distribution  at 
t  =  0.01  are  shown  in  Fig.  5.1.  The  macro  solver  from  left  to  right  in  Fig.  5.1 
is  the  DG  DG  P^  and  the  second-order  central  scheme,  respectively. 
We  can  see  that  the  DG  scheme  captures  the  shocks  much  better  than  the 
second-order  central  scheme  using  the  same  mesh.  We  can  also  see  this  from 
the  cross  section  view  of  the  3D  temperature  distribution  by  the  DG  and  the 
second-order  central  scheme  in  Fig.  5.2.  Since  we  have  not  used  any  limiter 
for  the  results  in  the  left  picture  of  Fig.  5.2,  there  are  some  oscillations  near 
the  shocks.  The  application  of  limiters  can  eliminate  these  oscillations,  see 
Fig.  5.2  right.  For  both  hgures,  the  shocks  are  sharper  for  the  DG  method 
than  for  the  second-order  central  scheme. 

With  the  heat  flux:  Adding  the  heat  flux  in  our  model  is  adding  a 
diffusion  term  to  the  right  side  of  PDFs.  The  LDG  methods  can  easily  deal 
with  the  diffusion  term.  The  results  of  the  temperature  distribution  using 
both  DG  P^  and  P^  at  t  =  0.01  with  a  80  x  80  mesh  are  shown  in  Fig.  5.3. 

The  full  MD-based  constitutive  model 
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Table  5.3:  error  and  order  of  accuracy  for  the  energy 

mesh  standard  local  curl-free 

error  order  L^error  order 


pi 

10  X  10 

3.38E-02 

- 

3.38E-02 

~ 

20  X  20 

7.97E-03 

2.08 

7.98E-03 

2.08 

40  X  40 

7.42E-04 

3.43 

7.42E-04 

3.43 

80  X  80 

1.75E-04 

2.08 

1.75E-04 

2.08 

160  X  160 

4.42E-05 

1.98 

4.42E-05 

1.98 

p2 

10  X  10 

1.34E-01 

- 

1.12E-01 

~ 

20  X  20 

1.78E-02 

2.91 

1.54E-02 

2.87 

40  X  40 

2.29E-03 

2.96 

1.99E-03 

2.95 

80  X  80 

2.92E-04 

2.97 

2.54E-04 

2.97 

160  X  160 

3.68E-05 

2.99 

3.21E-05 

2.98 

In  this  case,  we  compute  the  stress  from  MD  in  every  cell.  We  first 
compute  the  temperature  from  the  atomistic  models.  After  the  temperature 
is  obtained,  we  compute  the  stress  and  energy  flux  from  the  atomistic  model 
via  a  canonical  ensemble. 

Without  the  heat  flux.  The  results  of  the  temperature  distribution  at 
t  =  0.01  are  shown  in  Fig.  5.4.  We  compare  our  DG  method  again  with  the 
second  order  central  scheme.  Using  four  times  as  many  cells  for  the  second- 
order  central  scheme  as  that  for  the  DG  method,  we  can  still  see  the 
advantage  of  the  DG  method  in  capturing  shocks.  A  clearer  view  in  Fig.  5.5 
shows  the  cross  section  of  the  3D  temperature  distribution  computed  by  the 
DG  and  the  second-order  central  scheme.  The  results  of  the  velocity  held 
at  t  =  0.01  are  shown  in  Fig.  5.6.  The  direction  of  the  velocity  is  pointing 
outward  since  the  temperature  is  propagating  outward. 

With  the  heat  flux.  The  results  of  the  temperature  distribution  and 
the  velocity  held  at  t  =  0.01  are  shown  in  Fig.  5.7. 

The  domain  decomposition  model 

Here  we  use  the  DDM  procedure,  as  described  in  Sec.  5.2. 
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Without  the  heat  flux:  The  results  of  the  temperature  distribution 
and  the  velocity  flelds  using  DG  aX  t  =  0.01  by  the  DDM  model  are 
shown  in  Fig.  5.8. 

The  results  by  the  DDM  are  almost  the  same  as  the  one  we  obtained  by 
the  full  MD-based  model  (see  Fig.  5.4  left  and  5.6  left).  This  is  more  clearly 
seen  from  the  cross  section  of  the  3D  temperature  distribution  by  the  DDM 
and  the  MD  (see  Fig.  5.11):  they  overlap  almost  completely. 

The  results  of  the  temperature  distribution  and  the  velocity  flelds  using 
the  DG  at  t  =  0.03  by  the  DDM  model  are  shown  in  Fig.  5.9.  As  the 
time  evolves,  the  heating  is  expanding  outwards. 

With  the  heat  flux:  The  results  of  the  temperature  distribution  and 
the  velocity  flelds  at  t  =  0.01  are  shown  in  Fig.  5.10.  Fig.  5.12  shows  the 
results  at  t  =  0.03.  We  can  clearly  see  a  thermal  expansion  propagating 
outwards. 

The  domain  decomposition  method  (DDM)  saves  lots  of  computational 
cost.  For  example,  we  need  18  hours  using  36  processors  to  evolve  1  time 
step  with  30  x  30  cells  for  the  MD  DG  P^  method,  while  we  only  need  less 
than  3  hours  to  do  the  same  thing  for  the  DDM.  Actually,  we  only  need  to 
run  the  MD  for  the  high  temperature  region  in  the  middle.  Initially,  this 
part  is  only  16%  of  the  whole  domain.  As  time  evolves,  the  high  heat  will 
propagate  outwards  and  the  area  becomes  even  smaller. 

5.4.3  Wave  propagation 

As  the  last  example,  we  simulate  elastic  wave  propagation  in  a  2-D  un¬ 
bounded  medium.  The  experiment  we  present  is  described  in  [5].  The  major 
difference  is:  their  governing  equations  are  based  on  linear  elastodynamics, 
but  ours  are  based  on  the  domain  decomposition  model.  We  apply  perfectly 
matched  layers  (PML)  as  described  in  [18]  to  simulate  the  propagation  of 
waves  in  this  open  domain.  See  Appendix  B  for  a  brief  description  of  the 
PML  method. 

Our  computational  domain  is  D  =  [—5,  5]  x  [—5,  5]  in  2-D,  occupied  by 
an  elastic  material,  and  we  suppose  that  the  initial  condition  (or  the  source) 
is  supported  in  D  =  [—3.5,  3.5]  x  [—3.5,  3.5]  C  D.  We  are  solving  this  elas¬ 
todynamics  problem  with  absorbing  layers  (PML)  with  width  5  =  1.5  on  all 
four  boundaries. 

The  initial  data  of  the  velocity  and  the  stress  are  taken  to  be  equal  to 
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zero.  The  material  is  of  constant  temperature  T  =  0.167  everywhere,  and 
it  is  assumed  to  stay  in  the  same  temperature  environment.  An  explosive 
source  located  at  the  point  {xs,ys)  =  (—3.15,3.15)  is  introduced  to  the  right 
side  of  the  equations  of  conservation  of  velocity. 


=  h{t)g{r)er. 


(5.16) 


where  e).  =  (x  —  Xgiy  —  ys),  r  =  ||e).||.  The  function  h(t)  is  the  so-called 
second-order  Ricker  signal  with  central  frequency  equal  to  /o  =  lOHz  (see 
Fig.  5.13  left): 

h{t)  =  [2n\fot  -  1)2  -  (5.17) 

and  the  function  g{r)  is  the  Gaussian  function  dehned  by  (see  Fig.  5.13  right) 


9{r) 


10e-7hAo)2 


r 


2 

0 


1 


(5.18) 


which  is  concentrated  in  a  small  disk  of  radius  tq  =  0.3. 
The  damping  factor  we  use  is  as  follows: 


d{x) 


(5.19) 


where  R  =  10“^  is  the  theoretical  reflection  coefficient  from  the  terminating 
reflection  boundaries,  i.e.,  the  reflection  coefficient  is  about  0.1%.  c  is  an 
upper  bound  of  the  wave  velocities,  chosen  as  c  =  1  here. 

We  show  the  snapshots  of  the  velocity  module  using  DG  by  the  macro 
model  in  Fig.  5.14  and  by  the  DDM  model  in  Fig.  5.15.  The  PML  works 
very  well.  The  wave  fronts  are  not  circles  because  of  the  anisotropy  of  the 
medium.  The  waves  are  propagating  outwards.  At  t  =  0.34  (see  Fig.  5.14, 
5.15  upper  right),  some  waves  disappear  on  the  upper  left  corner.  The  lower 
right  waves  are  moving  downwards  as  the  time  goes.  At  t  =  1.0  (see  Fig.  5.14, 
5.15  bottom  middle),  these  waves  arrive  at  the  lower  right  corner.  And  at 
t  =  1.5  (see  Fig.  5.14,  5.15  bottom  right),  the  amplitude  of  velocity  module 
is  actually  already  below  5  x  10“"^. 


6  Concluding  remarks 

We  have  developed  a  multiscale  solver  based  on  the  discontinuous  Galerkin 
(DG)  method.  The  ability  of  the  DG  method  to  treat  sharp  wave  front  makes 
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it  particularly  suitable  for  this  problem.  This  has  been  demonstrated  by  a  few 
test  problems.  More  importantly  it  allows  us  to  conveniently  use  atomistic 
models  to  provide  constitutive  data  within  the  computation.  The  bottleneck 
of  the  method  has  been  the  atomistic  component  where  MD  is  performed  to 
calculate  the  elastic  stress.  This  has  been  alleviated  to  some  extent  in  this 
paper  by  hrst  obtaining  a  simplihed  model  in  the  case  of  small  strain  and  low 
energy/temperature.  More  efficient  methods  are  needed  here  to  completely 
bypass  MD.  This  is  our  current  work  in  progress.  Finally  even  though  the 
numerical  tests  have  been  performed  using  rectangular  elements,  the  DG 
methodology  allows  the  usage  of  arbitrary  unstructured  meshes  with  a  full 
h-p  adaptivity  capability. 

In  the  current  framework,  we  have  assumed  that  the  material  is  a  single 
perfect  crystal.  The  next  step  is  to  model  the  dynamics  of  isolated  de¬ 
fects,  such  as  dislocation,  phase  boundary  or  crack  tips.  These  problems 
bring  up  many  interesting  issues,  such  as  the  boundary  conditions  at  the 
atomistic/continuum  interface  (e.g.  see  [37,  38]),  adaptive  mesh  rehnement, 
thermal  fluctuations  etc.  These  issues  will  be  addressed  in  future  work. 


A  Appendix:  The  hyperbolicity  of  the  con¬ 
stitutive  model  without  the  heat  flux 

When  the  deformation  is  small  and  the  temperature  is  low,  we  can  use  the 
linearized  stress  strain  relation  (5.4)  and  the  linearized  temperature  relation 
(5.6)  to  solve  for  the  temperature  T,  then  we  will  have  a  closed  system  of 
equations.  For  the  case  without  the  heat  flux,  this  system  of  equations  can 
be  written  in  the  form  of  a  system  of  conservation  laws: 

Ut  -|-  Fj,  -|-  =  0, 

where 


(  ^ 

/  -vi  \ 

/  0 

^12 

0 

-Vi 

^21 

-V2 

0 

^22 

,  F  = 

0 

,  G  = 

-V2 

P0^'1 

-Pll 

-P2I 

P0V2 

-P21 

—P22 

\  Poe  +  \po{vf  +  v^)  y 

\  -ViPu  -  V2P21  ) 

\  -V1P2I  -  ' 

(A.l) 
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with  the  relations  (5.4)  and  (5.6). 

It  is  known  that  in  [20]  (p.55),  this  system  (A.l)  is  hyperbolic  on  a  certain 
region  of  the  state  space  if  for  every  (e,  77)  lying  in  that  region,  the  following 
two  conditions  hold, 


dE{e,ri) 

dr} 


>0, 


(A.2) 


>  0,  for  all  V  and  f,  (A.3) 

OEijOeki 

where  77  is  the  entropy  and  E{e,  if)  =  E{e,  T).  The  first  condition  ensures  that 
the  temperature  is  positive,  and  the  second  condition,  often  referred  to  as 
the  rank-one  convexity  condition,  implies  that  the  local  thermal  equilibrium 
is  stable. 

First  to  compute  the  entropy  77  =  r]{e,T),  using  the  Helmholtz  free  energy 
([20]  p.41) 

if{e,T)  =  e{e,T)-T7},  (A.4) 

we  obtain 


Eij  Po 


dfj  dE{e,  T) 


de 


de. 


dr}  rp 

PoT-e —  =  CijkiSki  —  PoT 


de. 


de 


(A.5) 


The  last  equality  comes  from  the  expression  of  E{e,T)  given  by  (5.6). 
Applying  the  dehnition  of  P  in  (5.3)  to  (A.5)  ,  we  get 

PoT-^  =  -aT5,,. 


de. 


Thus 


a 


V  = - (£11  +  ^22)  +  ^{T). 

Po 


(A.6) 

(A.7) 


To  determine  the  function  v?(T),  we  know  from  (A.4), 


,  dp 


rr\  _  dfj  _  d{e-Tr})  _  de  ^  ... 

pi  ^  dT^^^^ dT  ^  ^ 


which  then  reduces  to 


dT  Po' 


(A.9) 
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Together  with  (A. 7),  we  hnally  obtain  the  entropy  for  the  system  (A.l), 


r]{e,T)  =  -  —  {eii  +£22)  +  —  logT 
Po  Po 

up  to  a  constant. 

Next  ,  f  =  f{e,ri)  can  be  solved  by  (A. 10), 


(A.IO) 


Then 


=  g(P0)?+a(£ll+£22))/3_ 


For  the  first  condition  given  by  (A. 2), 


dE{e,  if)  _  ^^^(^p^r^j^a(eii+e22)) /i  ^  g_ 
dr] 


(A.ll) 

(A.12) 

(A.13) 


The  second  condition  (A. 3),  called  the  Legendre-Hadamard  condition, 
means  that  E  is  rank-one  convex  in  £.  It  is  equivalent  to  the  condition  that 

the  matrix  (  )  is  positive  dehnite. 

\deijdeki  J  n 

It  is  easy  to  check  that  the  matrix 


7  d^E{£,ri) 

\deijdeki  ) 

(  Cn  +  ia2g(po>?+aFii+£22))/3  Q  Q  C12  +  la2e(^0’?+“(^ii+^22))/3 

0  C*44  C44  0 

0  C*44  C44  0 

V  Cii  +  0  0  Cii  +  j 

(A.14) 

is  positive  dehnite. 

Therefore,  the  system  (A.l)  is  hyperbolic. 


B  Appendix:  PML  methods 

The  idea  of  the  PML  is  to  surround  the  computational  domain  with  an 
absorbing  layer  (the  PML  region)  such  that  the  coupled  system  possesses  the 
property  of  generating  no  rehection  at  the  interface  between  the  free  medium 
and  the  artihcial  absorbing  medium.  The  principle  ideas  of  PML  has  been 
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first  introduced  by  Berenger  [6].  In  [5],  the  authors  have  shown  that  it  is 
possible  to  write  in  a  systematic  way  of  designing  a  PML  model  for  a  general 
hrst-order  hyperbolic  system.  In  this  section,  we  will  briefly  describe  the 
general  construction. 

Consider  a  general  two  dimensional  first-order  hyperbolic  system 

Ut  +  AUx  +  BUy  =  0,  I'R  1 'i 

Uit  =  0)  =  Co,  ^  ^ 

where  C  is  a  m— vector,  A  and  B  are  m  x  m  matrices.  To  simplify  the 
presentation,  we  assume  our  physical  domain  is  in  the  left  half  plane  i.e. 
[xi  <  0,a;2  =  0]  X  [|/i,|/2].  The  initial  condition  Uq  is  zero  on  the  right  half 
plane.  The  main  idea  of  the  PML  model  is  to  couple  the  equation  in  the 
left  half-space  with  an  equation  in  the  right  half  space  such  that  there  is  no 
reflection  at  the  interface  y  =  0. 

The  construction  of  PML  in  the  x-direction  is  to  split  U  =  -|-  U^‘^\ 

such  that  the  unknown  IjA'^  is  only  associated  to  the  derivatives  with  respect 
to  x,  and  to  the  derivatives  with  respect  to  y.  Then  we  introduce  a 
damping  factor  d{x)  only  on  We  obtain  the  system  in  the  following: 

+  d{x)U^^^  +  AUx  =  0, 

ui^^+BUy  =  0,  (B.2) 

C(t  =  0)  =  Co, 

where  d{x)  =  0  for  a;  <  0,  d{x)  <  0,  for  a;  <  0.  It  is  easy  to  see  that  C 
satishes  the  same  system  of  equations  in  the  physical  domain.  The  analysis 
[6]  shows  that  there  will  be  no  reflection  at  the  interface  between  the  phys¬ 
ical  domain  and  the  absorbing  layers.  Furthermore,  the  transmitted  wave 
decreases  exponentially  during  its  propagation  inside  the  layer. 

In  our  example,  the  macro  equations  are  a  2D  nonlinear  system  of  parabolic 
equations.  We  will  apply  the  PML  to  all  the  boundaries.  So  our  macro  equa¬ 
tions  can  be  rewritten  in  the  following  form; 

+  d{x)U^^'^  +  Fx  =  0, 

Cf^ +d(l/)C(2) +Gy  =  0,  (B.3) 

C(f  =  0)  =  Co, 

where  F  and  G  are  in  (A.l),  d{x)  and  d{y)  are  zeros  in  the  physical  domain 
and  positive  in  the  absorbing  layers. 
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(a)  A  conventional  numerical  method. 


(b)  The  heterogeneous  multiscale  method. 

Figure  3.1:  Conventional  numerical  methods  and  the  heterogeneous  multi¬ 
scale  methods.  In  the  multiscale  method,  an  additional  component,  indicated 
by  the  small  boxes,  is  applied  to  compute  the  elastic  stress.  Here  9  =  ksT. 
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Figure  4.1:  Numerical  test  on  ID  shock  formation  and  propagation.  100 
macro  cells  are  used.  The  solutions  are  displayed  at  the  physical  time  T  = 
0.01.  Solid  line  with  the  square  symbol:  the  second-order  central  scheme; 
triangle  and  circle  symbols:  DG  and  with  TVD  limiter.  Left:  strain, 
right:  velocity. 
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Figure  5.1:  Temperature  distribution  without  the  heat  flux  (by  the  simplihed 
constitutive  model)  at  f  =  0.01  with  80  x  80  cells.  Left:  DG  P^  with  the 
TVD  limiter;  middle:  DG  P^  with  the  TVB  limiter  M  =  10;  right:  the 
second-order  central  scheme. 
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temperature 


Figure  5.2:  Cross  section  of  the  3D  temperature  distribution  without  the  heat 
flux  (by  the  simplihed  constitutive  model)  at  t  =  0.01  with  80  x  80  cells.  Left: 
DG  without  limiters  and  the  second-order  central  scheme;  Right:  DG 

with  the  TVD  limiter  and  F^  with  the  TVB  limiter  M  =  10. 


DG  P1  N=80  t=0.01  with  heat  flux  DG  P2  N=80  t=0.01  with  heat  flux 


Figure  5.3:  Temperature  distribution  with  the  heat  flux  (by  the  simplihed 
constitutive  model)  at  t  =  0.01  with  80  x  80  cells.  Left:  DG  F^;  right:  DG 
FT 
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DG  P1  N=30  m=0 


2nd  central  scheme  N=60 


Figure  5.4:  Temperature  distribution  without  the  heat  flux  (by  the  full  MD- 
based  model)  at  t  =  0.01.  Left:  DG  with  30  x  30  cells  with  the  TVD 
limiter;  right:  the  second-order  central  scheme  with  60  x  60  cells. 


Figure  5.5:  Cross  section  of  the  3D  temperature  distribution  without  the 
heat  flux  (by  the  full  MD-based  model)  at  f  =  0.01.  DG  with  30  x  30 
cells  and  the  second-order  central  scheme  with  60  x  60  cells. 
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Figure  5.6:  The  velocity  field  without  the  heat  fiux  (by  the  full  MD-based 
model)  at  t  =  0.01.  Left:  DG  with  30  x  30  cells  with  the  TVD  limiter; 
right:  the  second-order  central  scheme  with  60  x  60  cells. 


DG  P1  N=30  T=0.01  with  heat  flux 


DG  P1  N=30  T=0.01  with  heat  flux 


Figure  5.7:  With  the  heat  fiux  at  t  =  0.01.  DG  with  30  x  30  cells  (by 
the  full  MD-based  model).  Left:  the  temperature  distribution;  right:  the 
velocity  field. 
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Figure  5.8:  Without  the  heat  flux  at  t  =  0.01.  DG  with  30  x  30  cells  (by 
the  DDM).  Left:  the  temperature  distribution;  right:  the  velocity  held. 
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Figure  5.9:  Without  the  heat  flux  at  t  =  0.03.  DG  with  30  x  30  cells  (by 
the  DDM).  Left:  the  temperature  distribution;  right:  the  velocity  held. 
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DDM  DG  P1  N=30  T=0.01  with  heat  flux 


Figure  5.10:  With  the  heat  flux  at  t  =  0.01.  DG  with  30  x  30  cells  (by 
the  DDM).  Left:  the  temperature  distribution;  right:  the  velocity  held. 


Figure  5.11:  Comparison  of  DDM  and  full  MD-based  models:  cross  section 
of  the  3D  temperature  distribution  at  t  =  0.01.  DG  with  30  x  30  cells. 
Left:  without  the  heat  flux;  Right:  with  the  heat  flux. 
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DDM  DGP1  N=30  T=0.03  with  heat  flux 


DDM  DGP1  N=30  T=0.03  with  heat  flux 
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Figure  5.12:  With  the  heat  flux  at  t  =  0.03.  DG  with  30  x  30  cells  (by 
the  DDM).  Left:  the  temperature  distribution;  right:  the  velocity  held. 


Figure  5.13:  Left:  the  Ricker  signal  h(f);  right:  the  Gaussian  function  g{r). 
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T=0.25 


T=0.34 


Figure  5.14:  Snapshots  of  the  velocity  module  at  different  times  (by  the 
simplihed  constitutive  model):  DG  with  80  x  80  cells. 
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Figure  5.15:  Snapshots  of  the  velocity  module  at  different  times  (by  the  DDM 
model);  DG  with  80  x  80  cells. 
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